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This work illustrates the optimization of thermoacoustic systems, while taking into account thermal 
losses to the surroundings that are typically disregarded. A simple thermoacoustic engine is used as an 
example for the methodology. Its driving component, the thermoacoustic regenerator (also referred to as 
the stack), is modeled with a finite element method and its dimensions are varied to find an optimal 
design with regard to thermal losses. Thermoacoustic phenomena are included by considering acoustic 
power, and viscous and capacitive losses that are characteristic for the regenerator. The optimization 
considers four weighted objectives and is conducted with the Nelder-Mead Simplex method. When 
trying to minimize thermal losses, the presented results show that the regenerator should be designed to 
be as short as possible. It was found that there is an optimal regenerator diameter for a given length. The 
results are presented for a variety of materials and weights for each objective. 

© 2009 Elsevier Masson SAS. All rights reserved. 


1. Introduction 

This work shows how the principles of optimization can aid in 
the design process of thermoacoustic devices. Thermoacoustic 
devices utilize sound waves to drive a thermodynamic process 
instead of mechanical pistons. Their advantage is the inherent 
mechanical simplicity. While the concept is not new, the tech¬ 
nology has not been advanced to a high degree, as compared to, for 
example, the internal combustion engine. Specifically, optimization 
has not been utilized well as a design tool. Rather, designs of 
thermoacoustic devices rely on parametric studies and “rule-of- 
thumb” experience. In this context, optimization will be defined as 
modeling the entire system and allowing all degrees of freedom to 
be varied as part of the investigation, as opposed to the parametric 
study, where only one parameter is varied while all others are kept 
constant. Before elaborating of the optimization of a thermoacous¬ 
tic device, we will first introduce the concepts of thermoacoustics 
in detail. This will also illustrate some of the terminology used in 
the sections discussing the model, the objective functions and the 
results. 

The basic thermodynamic cycle occurring in thermoacoustic 
devices is the Stirling cycle, which was developed in 1816 [1]. The 
original mechanical Stirling engine utilized two pistons and 
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a regenerative heat exchanger [2]. Over the course of one cycle, 
the working gas is compressed, and it then transfers thermal 
energy to the heat sink, thus maintaining a constant temperature. 
Afterwards, the gas is heated at constant volume by the regen¬ 
erator and then is heated further at the heat source. This heat 
supply occurs while the gas is allowed to expand and drive the 
power piston, again at constant temperature. After expansion, the 
gas is displaced to the heat sink, while cooling off at constant 
volume by depositing heat to the regenerator, which stores heat 
between cycle segments [2]. It is noteworthy that this externally 
heated, closed cycle uses the same gas for all stages, as opposed 
to the internal combustion engine, which has a constant 
throughput of working gas and fuel. The first application of this 
cycle as a thermoacoustic technology occurred when Ceperley 
recognized that sound waves could replace pistons for gas 
compression and displacement [3]. Since then, a wide variety of 
thermoacoustic engines (TAEs) and their counterpart, thermoa¬ 
coustic refrigerators (TARs), have been developed. Engines utilize 
a heat input to create intense sound output, while refrigerators 
can utilize this intense sound to withdraw energy form their 
surroundings. 

1.2. Thermoacoustic engines 

There are two main approaches to thermoacoustic engines, 
namely standing wave and traveling wave devices. Both contain 
a regenerative unit (called a stack or regenerator, respectively) 
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Nomenclature 

lx) 

angular frequency (s -1 ) 


area (m 2 ) 

n 

perimeter (m) 

A 

VT 

temperature gradient (Km -1 ) 

c 

speed of sound (ms -1 ) 



C 

capacitance (m -1 ) 

Dimensionless groups 

c p 

heat capacity (J kg -1 K -1 ) 

Gr 

Grashof number 

f 

frequency (s -1 ) 

Nu 

Nusselt number 

g 

gravitational acceleration 

Pr 

Prandtl number 

h 

heat transfer coefficient (W m -2 1< -1 ) 

Ra 

Rayleigh number 

H 

height (m) 



k b 

Boltzmann constant 

Subscripts and superscripts 

k 

thermal conductivity (W m -1 1< -1 ) 

1 

first order 

l 

plate thickness (m) 

00 

ambient, free stream 

L 

inertance (kg m -4 ), length (m) 

C 

channel 

N 

number of 

char 

characteristic 

P 

pressure (N m -2 ) 

crit 

critical 

Q 

heat flux (W) 

cold 

cold side 

R 

resistance (kg m -2 s -1 ) 

cond 

conductive 

T 

temperature (K, °C) 

conv 

convective 

u 

velocity (ms -1 ) 

D 

diameter 

W 

acoustic work (W) 

hot 

hot side 

y 

plate spacing (m) 

K 

thermal 



m 

time averaged 

Greek symbols 

obj 

objective 

d 

penetration depth (m) 

rad 

radiative 

£ 

surface emissivity coefficient 

s 

solid, standing 

£ 

plate heat capacity ratio 

xx,xy,yx,yy tensor directions 

T 

isentropic coefficient 

V 

viscous 

r 

temperature gradient ratio 

w 

wall 

p 

density (kg m -3 ) 




sandwiched between two heat exchangers, one to supply heat at 
high temperature (on the order of several hundred degrees Celsius), 
the other to withdraw heat from the system at (ideally) ambient 
temperature. In practice, the cold side has to be cooled because of 
conduction of heat from the hot side to the cold side, thus heating 
this side to temperatures higher than ambient. The temperature 
gradient across this porous section results in amplification of 
pressure disturbances in the working gas and results in a loud noise 
being emitted once a steady state has been achieved. In order for 
amplification to occur, this temperature gradient has to be larger 
than the critical temperature gradient, which is related to the 
temperature gradient that the gas would experience if it were 
under the influence of a sound wave in adiabatic conditions. The 
expression for this critical temperature gradient was derived by 
Swift [4] and is given in Equation (1.1): 



_wpi_ 

p m Cpll] 


( 1 . 1 ) 


This critical temperature gradient depends on the operating 
frequency w, the first order pressure and velocity in the standing 
wave pi and u\, as well as the mean gas density p m and specific heat 
c p . Depending on the ratio of the temperature gradient and critical 
temperature gradient, acoustic work is either created (a thermoa¬ 
coustic engine) when dr/dx* rit > 1 or transformed into heat energy 
(a refrigerator) when dT / dx < \ [5]. Fig. 1 shows a demonstration 
engine. It is a standing wave, quarter-wavelength device. The 
porous stack is located near the closed end of the resonance tube 
which is the pressure antinode and the velocity node. As a result, 
the gas inside the stack experiences large pressure oscillations and 


relatively small displacement. The heat input is provided with 
a heating wire. Active cooling with a second heat exchanger is not 
necessary because of the low thermal conductivity of the ceramic 
regenerative unit. 

The phasing of a standing wave is such that pressure and 
velocity are out of phase, and as a consequence, the heat transfer 
between the gas and the wall has to be delayed artificially. This 
delay is achieved by utilizing the poor thermal contact between gas 
and solid. As a result, however, this heat transfer also causes 
significant entropy generation, which limits efficiency. In order to 
avoid these losses, practical device designs have utilized traveling 
wave phasing, where velocity and pressure are in phase. In this 
case, the thermal contact between gas and solid does not have to be 
delayed; the efficiency of such traveling wave devices is inherently 
better than that of standing wave devices. This traveling wave 
phasing between pressure and velocity can be achieved in several 
ways. For example, one can use a looped tube as the compliance or 



Fig. 1 . Picture of a simple standing wave engine demonstrator. 
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one can utilize an annular so-called feedback inertance. The looped 
tube is illustrated by Swift in [6]. The latter has been utilized by 
Bastyr et al. (refer to [7] and [8]) and is very similar to a standing 
wave device in appearance. Both create a shift in the phase of the 
pressure wave relative to the velocity wave so that they are actually 
in phase with each other. This allows for the maximum displace¬ 
ment to occur at the moment of greatest compression or rarefac¬ 
tion, which is advantageous to the underlying thermodynamic 
cycle. Backhaus et al. introduced a large scale traveling wave engine 
that achieved a thermal efficiencies of 30%, corresponding to 41% of 
Carnot efficiency [9]. In any device, the resonator determines the 
operating frequency [10]. 

To illustrate the thermoacoustic phenomenon, we can draw 
parallels to optics. The amplification of the acoustic wave is 
similar to an optical laser, where light waves travel between 
a mirror and a partially silvered mirror in a standing wave 
fashion. The light waves are amplified through resonance and 
released through the partially mirrored side as a high power laser 
beam. Similarly, the amplified sound waves can also be extracted 
from the resonance tube to power external devices [11]. 
Commonly, the only application of thermoacoustic engines is to 
drive TARs, which utilize the reverse Stirling cycle, that attenuate 
the pressure in a sound wave to withdraw heat from the 
surroundings. 

2.2. Thermo acoustic refrigerators without moving parts 

Refrigerators and chillers driven by TAEs are a reality today; 
however they are limited to few, specialized uses, such as gas 
liquefaction. There are several advanced models of TARs used for 
these specialized applications. Radebaugh at the National Insti¬ 
tute of Standards and Technology (NIST) built a TAR with 5 W of 
cooling power at 120 I< and a low temperature at no load of 90 I< 
[12]. The main benefit is the lack of moving parts (such as seals), 
thus reducing maintenance costs. It is noteworthy that TARs can 
achieve these low temperatures in a single stage, whereas 
conventional vapor compression refrigerators (VCRs) can only 
achieve approximately 2301< in a single stage [12]. Cryogenic 
cooling is not the only application for TARs. A practical example 
of this design was given by Poese with a freezer for ice cream 
storage. This small scale chiller featured an annular space around 
the regenerative unit to achieve traveling wave phasing [8]. To 
date, the applications of thermoacoustics to the field of refriger¬ 
ation has been very limited, while the advantages have been 
demonstrated. In addition to the mechanical advantages, TARs do 
not use harmful refrigerants and can be driven by waste heat. The 
former is important with respect to the environmental concerns 
of refrigeration in general, while the latter is an important aspect 
when refrigeration has to be provided in locations where elec¬ 
trical power is not available or could be conserved. Consider, for 
example, the air-conditioning system in a passenger vehicle. The 
engine gives off large amounts of waste heat at high temperature, 
which could be utilized by a TAR instead of a mechanical air- 
conditioning system. The result is reduced fuel consumption 
because the TAR based system does not require pumps or 
compressors. 

2.3. Previous optimization efforts 

Currently, thermoacoustic technology is not as advanced as the 
internal combustion engine, which has experienced significant 
improvement since its conception over a century ago. As a result, 
there is much room for discovery in the field of thermoacoustics. 
Optimization techniques as a design aid, for example, are severely 
under-utilized, and previous efforts in thermoacoustic optimization 


are rare. Minner et al. consider the optimization of a thermoa¬ 
coustic refrigeration system. This work uses extensive model 
development and seeks to optimize the coefficient of performance. 
The group considers geometric parameters and fluid properties of 
the system and a simplex algorithm to search for the optimal 
solution. However, in order to account for the thermoacoustic 
operating conditions, DeltaE is used extensively [13]. Both Wetzel 
and Besnoin discuss optimization of thermoacoustic devices in 
their work. Wetzel targets the optimal performance of a thermoa¬ 
coustic refrigerator, Besnoin targets the heat exchangers [14] and 
[15], respectively. 

In addition to these optimization efforts, parametric studies 
have been utilized to estimate the effect of single design 
parameters on device performance. Zoontjens illustrated the 
optimization of inertance sections of a thermoacoustic devices. 
Upon closer inspection, they used DeltaE to vary individual 
parameters to determine optimal designs [16]. Ueda also deter¬ 
mined the effect of a variation of certain engine parameters on 
pressure amplitudes [17]. Tijani et al. attempted to optimize the 
stack spacing; however, they also utilized DeltaE for this work 
[18]. This is by no means a complete list of the “optimization” of 
engine components, but it is a good overview of optimization 
targets. Each work is undoubtedly a valuable addition to the 
thermoacoustic community, but they should not be considered 
optimizations in the classical sense, but rather parametric studies. 
In all likelihood, each “optimal” design is a local optimum as the 
optimization performed by each group considers one variable and 
all else equal. 

One common trait of all previous optimization efforts is that 
thermal losses to the surroundings that occur in the operation of 
the devices are not considered. In reality, however, we can optimize 
a thermoacoustic device with regard to one of the following 
objectives: 

1. Power output, 

2. Heat input, 

3. Viscous losses in the individual channels, and 

4. Heat loss through the device boundaries and cooling medium. 

It is intuitive that an increased heat input may increase the 
power output, but as a result, it may also increase the heat losses. It 
is the goal of this work to illustrate these tradeoffs, as they are not 
understood in detail. Below, an estimate of the magnitude of these 
thermal losses will be provided. It shows that these losses are 
significant when compared to total heat input and should be 
considered as a design criterium. 

2.4. Motivation 

Considering a simple thermoacoustic engine, comprised of 
a stack inside of a resonance tube, the energy flows are obvious. 
One side of the stack is heated, while the other side is cooled. As 
a result of this heating, strong pressure oscillations are created, 
which is the expected output of any thermoacoustic engine. 
However, as a secondary result of this heating, we must consider 
the thermal losses in the stack region of the engine. All three heat 
transfer modes are present: convection, radiation and conduction. 
Below, an estimate of the magnitude of said energy losses is per¬ 
formed. The values of each parameter discussed are derived from 
measurements from a simple demonstration engine as shown 
above in Fig. 1. During these measurements, the electrical power 
supplied to the ceramic stack was increased until oscillations 
occurred. At the onset, it was recorded that the hot side tempera¬ 
ture was 306 °C and the cold side (uncooled, facing ambient air) 
showed a temperature of 50 °C. The applied voltage was 15.1 V. 
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With a known resistance of the heating wire, an input power of 
approximately 50 W can be calculated. Given an outer diameter of 
the glass tube that houses the stack, and a uniform temperature 
distribution across the entire cross section of the stack and tube, 
and a total length of 20 mm over which this temperature difference 
is applied, we can calculate a convective heat loss (free convection, 
using a temperature dependent Nusselt number) of 4.7 W. In 
addition, we have to consider a radiative heat loss as well. Using the 
surface coefficient of emission provided in the pyrex data sheet 
(^Pyrex — 0.85...0.95), this radiative heat loss can be estimated to be 
approximately 5.1 W. This estimate shows that the heat losses, for 
the case of the demonstrator engine, are significant relative to the 
total energy supplied to it. This provides the motivation to optimize 
the geometry of the stack. There are obvious tradeoffs between the 
length of the stack and the circumference of the tube that houses 
the stack. It is intuitive that an increased heat input may increase 
the power output, but as a result, it may also increase the heat 
losses. 

The magnitude of the thermal losses leads us to consider them 
as a design criterium for thermoacoustic devices. As a starting point 
this work target a thermoacoustic engine. These systems provide 
a simple geometry to demonstrate the feasibility of our method¬ 
ology. TAEs will ultimately be applied in TAR systems to achieve 
refrigeration with no moving parts, which will undoubtedly 
increase system complexity and the modeling efforts. For future 
thermoacoustic systems designs, the thermal losses should be 
taken into account, especially with regard to the miniaturization of 
thermoacoustic devices. As the overall size shrinks, the surface area 
to (active) volume of a thermoacoustic engine, for example, 
increases, which in turn leads to higher thermal losses. This work 
aims to highlight one methodology to incorporate thermal losses in 
a design process by combining the fields of thermoacoustics and 
optimization. 

2. Model development 

The challenge in this work was to account for the temperature 
distribution throughout the stack in an accurate way. The solution 
to a three dimensional heat conduction problem subject to 
convective and constant temperature boundary conditions requires 
significant efforts [19]. For our initial model, only the stack geom¬ 
etry is considered; the model does not consider any variation in 
operating condition or the interdependency of stack location and 
performance. 

2.2. Computational domain 

Because of the symmetry present in the stack, the problem 
can be reduced to a two dimensional domain, with two constant 
temperature boundaries, one convective boundary, and finally 


Convection/ 

Radiation 



Adiabatic 


Fig. 2. Illustration of computational domain and implemented boundary conditions. 


one adiabatic boundary, as shown in Fig. 2. In order to maintain 
account for thermoacoustic phenomena, the work flow and 
viscous resistance occurring in the stack are account for, albeit in 
basic form. 

The temperature distribution throughout the rectangular 
domain was calculated using COMSOL Multiphysics, a finite 
element solver. The simple geometry allowed for straightforward 
computations when building and meshing the domain as well as 
solving for the temperature distribution. In order to consider all 
relevant energy fluxes correctly, the resulting outer surface area as 
well as the face areas have to be calculated by integrating over the 
angular component (defined as cp) of the cylindrical coordinates. In 
addition to calculating the temperature distribution in the stack 
efficiently, COMSOL also allows for integration into Matlab, which is 
important for the ensuing optimization, where the domain size was 
varied repeatedly. 

The rectangular representation of the stack was meshed using 
triangular cells. As the domain size will be varied over the course of 
the optimization, we have to maintain equal cell density for all 
domain sizes. Thus, the cell size was chosen to be constant for all 
domains used (i.e. fixed cell count per unit area). The cell type was 
chosen to be Lagrange quintic in order to minimize numerical 
errors. Lower order cell types caused some instability in the 
temperature solution throughout the domain as a result of the 
transition from initialized distribution to the physical distribution 
as the calculation progressed. 

2.2. Boundary conditions 

The boundary conditions on the modeled rectangle are 
prescribed as follows (left vertical wall as number 1, going counter¬ 
clockwise): 

1. Constant temperature (Thot), 

2. Adiabatic boundary, modeling the axis of the cylindrical stack, 

3. Constant wall temperature (T co id), and 

4. Free convection and radiation to surroundings (at Too), with 
a temperature dependent heat transfer coefficient. 

The material properties of the solid are assumed to be constant. The 
temperature of the domain is initialized to 300 K. The channels of 
the stack are not modeled explicitly, but they are accounted for by 
assuming an anisotropic thermal conductivity. 

2.3. Anisotropic thermal conductivity 

In order to account for the anisotropic nature of the stack as 
a result of the channels, the model was given two different thermal 
conductivities, one for the axial direction (open channels and solid 
in parallel), and one for the radial direction (open channels and 
solid in alternating series). COMSOL allows for anisotropic material 
properties in the form of 

k =('k X 1^) (2- 1 ) 

V K y x K yy J 

With this option, we can derive an advanced description of the 
stack without actually modeling the channels individually. As 
a basis, we consider the bulk thermal conductivity of the solid 
material being used, in addition to the channel size with the gas’ 
thermal conductivity for the net cross-channel thermal conduc¬ 
tivity. Only the values for axial and transverse thermal conductivity 
k xx and k yy are of interest. They can be calculated as shown in 
Equations (2.2) and (2.3): 
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Fig. 3. Temperature differences between isotropic and anisotropic thermal conduc¬ 
tivity for Copper and PMMA, Lx H = 0.4 x 0.2. 


kxx = 


fw^solid dck 


c'vgas 


tw d.Q 


( 2 . 2 ) 
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■yy 


kso\idkgas(tw + d c ) 
^solid^c + kg as tw 


(2.3) 


In order to show the influence of the anisotropic material proper¬ 
ties, we evaluated the temperature field in the rectangular domain 
for two different materials (copper and polymethyl methacrylate 
(PMMA)), in order to show two different extremes of bulk thermal 
conductivity. Fig. 3 shows the results. 


3. Objective functions 


With the computational domain defined, the objective func¬ 
tions can now be developed. In the present work, these are the 
conductive heat flux from the stack’s outer surface, the 
conduction through the stack as well as acoustic work and 
viscous resistance. 


3.1. Convective and radiative heat flux 

In the introduction, we have illustrated the magnitude of the 
convective heat flux from the outer surface of the stack to the 
surroundings in order to provide a motivation for this work. 
There, a linear temperature profile was used to estimate the 
heat transfer coefficient and the heat flux to the surroundings. 
In this model, however, the actual temperature distribution 
throughout the stack is be taken into account by utilizing 
COMSOL, as we account for a temperature dependence of the 
heat transfer coefficient. For the solid portion of the stack, these 
considerations assume bulk material properties that are inde¬ 
pendent of temperature. The convective heat transfer coefficient 
and the radiative heat flux to the surroundings are assumed to 
be dependent on the temperature. The total convective heat 
transfer across the cylindrical shell in its integral form can be 
described by: 


The heat transfer coefficient h is derived from the Nusselt 
number, which is a non-dimensional heat transfer coefficient. It is 
shown for the case of a horizontal tube subject to free convection 
[20]: 


Nu = 0.36 + 


0.518Ra^ 

~4 


i + ( 6 


kgas 


(3.2) 


This expression depends on the Prandtl number, a characteristic 
dimension D c har, which is the lateral domain dimension H in this 
case, and the Rayleigh number, which in turn can be expressed 
by: 


Ra 


Cr pr gg (T - Too )H 3 

VOL 


(3.3) 


where Gr is the Grashof number, Pr is the Prandtl number, T is the 
surface temperature, T x = 300 I< is the (constant) temperature of 
the surroundings, v is the viscosity of the surrounding gas, and a is 
the thermal diffusivity of the surrounding gas (air). Like all Nusselt 
number correlations, it is an empirical expression for this specific 
case of heat transfer. In order to derive the actual convective heat 
transfer coefficient, the Nusselt number is multiplied by the 
thermal conductivity of the surrounding gas, and divided by the 
characteristic dimension of the surface under consideration (again, 
stack radius H in this case). 

The radiation heat flux becomes increasingly important as Th 0 t 
increases, as shown in Equation (3.4): 


27r L 

Qrad = H/<b J f e(T(x) 4 -Ti)dxd<p (3.4) 

0 0 

The final heat flux objective for the top surface of the domain is the 
sum of both convective and radiative heat fluxes 


Qobj,l — Qconv + Qrad (3.5) 

where is the Stefan Boltzmann constant, and e is the surface 
emissivity, which depends on the emitted wavelength, and in 
turn is a function of temperature (an effect that may not be 
negligible if the temperature difference across the domain is 
sufficiently large). 


3.2. Conductive heat flux 

This heat flux is representative of the heat loss across the cold 
end of the domain. As the temperature gradient there is non-zero, 
a heat flux must be present. It is assumed that thermal energy is 
removed via the cooling water flow. Similar to the cylindrical shell, 
this heat flux has to be integrated over the whole surface repre¬ 
senting the cold side: 

2tt H 

Qcond = 0obj,2 = f f dX A<P (3-6) 

0 0 


2tt L 


Qconv = H JJ h(T(x))(T(x) -T^dxdtp. 


0 0 


In addition to the thermal energy fluxes inside the stack, we also 
have to account for energy fluxes that are inherent to thermoa¬ 
coustic engines. The acoustic power and the viscous resistance 
expressions are developed below. 
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3.3. Acoustic power 


According to Swift, the acoustic power per channel is given by 


W 


1 /477Lw 



(t - 1)p 2 

pc 2 ( 1 + e) 


(T - 1 ) - 5j,pu 2 


(3.7) 


which is composed of the thermal contribution minus the viscous 
effects [4]. 77 is the channel circumference, w is the operating 
frequency, y is the ratio of the specific heats (isochoric and 
isobaric), p is the gas’ density, and 7 is the ratio of critical 
temperature gradient and local temperature gradient, which has 
previously been introduced with Equation (1.1) [6]. p and u are the 
pressure and velocity amplitude achieved in the stack. With data 
taken from the engine assembly shown above, the pressure 
amplitude in the stack is assumed to be p max = 500 Pa. With this 
information, the velocity amplitude can be derived from Equation 
(3.8), which describes the velocity distribution in a standing wave: 


u = i^) cos j (3 - 8) 

to be Umax — Pmaxl(pc) — 1.3 ms -1 . For this study, both amplitudes 
are assumed to be constant for all stack geometries and materials. 
Finally, e is the so-called plate heat capacity ratio [21] and can be 
expressed by: 

_ (P C P^ K ) gas tan/i((i + l)yoAt) ,, 

0°^) S oiid tanft((i + l)l/a s ) 1 J 

This expression can be simplified to values of e=yold K if yold K < 1 
and e = 1 if yol5 K > 1, where yo is half of the channel height, / is half 
of the wall thickness, <5 S is the solid’s thermal penetration depth. As 
mentioned above, the expression for the work output is provided 
for a single channel. In order to provide a physical representation 
and to be consistent with the remaining assumptions, we estimate 
the total number of round channels N c as a function of cross section 
size for the cylindrical stack. For a circular domain, this estimate 
equals the ratio of overall domain size and cross sectional size of 
each channel, decreased to 70% of this value to account for the solid 
percentage of each channel. 



(3.10) 


An added degree of accuracy would be achieved by determining 
the actual number of channels by utilizing a packing correlation for 
square surfaces in a circular domain. For the sake of simplicity, this 
is omitted. 


3.4. Viscous resistance 


The viscous resistance for each channel is given by: 



pIIL 


(3.11) 


where 77 is the circumference and A c the area of the channel. This 
expression has the units [kg/m 4 s]. In order to express this in terms 
similar to the other variables used, we multiply Equation (3.11) by 
the volumetric velocity [m 3 /s] and the oscillating frequency [1/s], 
yielding [W/m] as a final unit for the viscous resistance per channel. 
Just as the total acoustic power of the stack was dependent on the 
total number of channels, the viscous resistance also depends on 
this value. As the full stack represents a network of parallel resis¬ 
tances, we divide the value for the individual resistance by the 


same factor N c derived above for the acoustic power. Below, we 
describe how the objective functions are conditioned and imple¬ 
mented in the optimization routine. 

4. Solution strategy 

Since we are attempting to find a geometric optimum of the 
stack geometry, the radius and length of the stack have to be varied 
as part of an optimization loop. Each time, the temperature distri¬ 
bution has to be calculated. The overall process can be summarized 
as follows: 

1. Initial guess for domain dimensions, 

2. COMSOL solution for temperature distribution and heat fluxes, 

3. Matlab evaluation of the objective and penalty functions, and 

4. Repeat if not optimal, else terminate. 


4.2. Normalizing objective functions 

The objective functions as introduced above in Section 3 are 
given in a general form, all with different units and also different 
orders of magnitude. This would result in skewed results from the 
optimization as the optimizer considers all objectives as a sum. For 
this reason, all objectives have to be scaled as to provide dimen¬ 
sionless values that are on a similar scale. 


4.2.2. Normalizing the work output and viscous resistance 

The maximum work output is always achieved when the 
domain size is maximized. This is obvious, both because the larger 
the cross section of the stack, the more channels it can hace, and 
because there is a linear relationship between power and stack 
length. Thus, the maximum work output is achieved at the 
maximum allowable length and cross section, multiplied by the 
open area coefficient, as defined in Equation (3.10), as shown in 
Equation (4.1): 


Wmax = N. 


cl H n 


d v 8 K (r — l)wp m ax f VT 
; pc 2 ( 1 + f ) \yr^i 


i 


(4.1) 


All parameters used here are the same as introduced above in 
Equation (3.7). For the viscous resistance, the maximum value is 
achieved with the longest channel, as this provides the most 
surface area and thus the highest viscous loss (Equation (4.2)). 
The maximum total resistance is achieved at the smallest tube 
cross section, because it can incorporate the smaller number of 
channels (which is accounted for by the factor N c | H . , in Equation 
(4.3)). 
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^,max — V^clH min ^,max,c 


(4.3) 


The work output and resistance expressions thus depend only on 
terms that are constant, and simple terms that depend solely on 
domain dimensions. For this reason, the normalization of both 
objectives is straightforward: 



value 

max.value 


(4.4) 


4.2.2. Normalizing the heat flux 

In the case of the heat fluxes, the normalization is less 
straightforward, as it is not immediately clear where their 
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maximum and minimum values are achieved. For this reason, 
a separate COMSOL calculation is performed, evaluating all heat 
fluxes for all possible combinations of domain dimensions that are 
allowed within the imposed boundaries on the domain (L m i n , L max , 
tfmin> tfmax)- This domain is discretized with a separate grid, and all 
heat flux values are calculated and stored in a matrix. The respec¬ 
tive maximum and minimum values of this matrix are then used to 
normalize the heat fluxes according to Equation (4.5): 


sufficiently large. On the other hand, a penalty that is too large 
can lead to problems when trying to solve the optimization 
problem; the problem can become ill-conditioned [22]. Conse¬ 
quently, the penalty function has to assume a value of 0 when 
the design constraints are not violated, and large if a bound is 
violated. The following sections will elaborate on the derivation 
of the objective functions and the penalty functions specific to 
this problem. 



value - (min.value) 
(max.value) - (min.value) 


(4.5) 


43. The optimization method: Nelder-Mead simplex 


With these normalization schemes in place, all objectives vary 
between 0 and 1, and both values are achievable within the 
domain. This allows us to optimize with respect to their sum, as 
the system is dimensionless and all four objectives are now on the 
same order of magnitude. In order to provide bounds on the 
variables used in this study, we implemented the following 
penalty functions. 

4.2. Penalty functions 

Two penalty functions are used for each decision variable (for 
the upper and lower bounds). For example, we want to avoid 
solutions such as the long and thin stack (“pencil” shape) or the 
short and large-radius stack (“disc” shape). The penalty functions 
are only evaluated if the variable is actually violating the boundary, 
otherwise they assume a value of 0. 

In our current case, the variables that need to be constrained are 
the axial and radial dimensions of the rectangular domain, L and H, 
respectively. They take the same format for both variables: 

if H_min - H > le-12 

penalty = 10 + (H_min - H) * le8 

end 

and 

if H - H_max > le-12 

penalty = 10 + (H - H_max) * le8 

end 

As soon as the constraint is violated, the constant term of the 
penalty ensures an immediate effect on the border of the 
feasible region, while the linearly increasing term ensures that 
the algorithm will avoid solutions further away from the feasible 
region. H m i n -H> le-12 rather than H m [ n -H>0 is used as 
a check for violation, as 0 is not well defined in numerical 
terms. The latter check could lead to numerical instabilities. In 
order to still provide a reasonable gradient as the domain 
violation becomes more severe, the difference H m [ n -H , and 
similarly for a violation of H-H max is multiplied by a large 
constant. The summation of the normalized objectives and the 
penalty functions form the total objective function F 0 bj as illus¬ 
trated in Equation (4.6). 

minimize F obj = w/ obj f + p, (4.6) 

i i 

where p z represent all penalty functions, and w/ is used to weight 
the individual objectives. The weighting factors can be used to shift 
the design emphasis on, say, the convective heat flux, or on work 
output, for example. 

It can be shown that unconstrained optimization using 
penalty functions can yield the optimal solution of the original 
objective if the penalty imposed on the new auxiliary function 
(which is the combination of primal function and the penalty) is 


The Nelder-Mead search algorithm in a d-dimensional space 
uses d +1 points to determine a downhill direction of an 
objective function. It does not rely on gradients, and thus the 
function does not have to be differentiable for the algorithm to 
be successful. For a surface, the algorithm uses three points 
(corners of a triangle), evaluates the function values and then 
applies one of four moves to the worst point: reflection, 
expansion, contraction, and shrinkage. The simplex changes its 
shape, and the one new function value is evaluated and 
compared to the one it has replaced. Based on the behavior of 
the function value from one iteration to the next, the algorithm 
chooses which modification to apply [23]. The algorithm uses an 
initial guess as input, and then varies the variables according the 
simplex function values. In the current case, it will vary both 
axial and lateral dimensions of the used rectangular domain. The 
hot side temperature is estimated before the first iteration based 
on the axial domain dimension and a given bulk temperature 
gradient. Locally, this temperature gradient varies significantly as 
the heat flux across the top surface is included (as calculated by 
COMSOL). 

5. Implemented optimization routine 

As introduced above, the optimization requires that we know 
the maximum possible value of each of the four objectives. For this 
reason, a preprocessor is implemented that utilizes COMSOL to 
calculate the heat fluxes across the top and cold surfaces for all 
possible domain permutations. The feasible reason is subdivided 
into a grid of L,H points and at each node, the objective functions 
are evaluated. The result is a matrix of objective function values as 
a function of both L and H. 

given Lmax, Hmax 

i j = 1 

for i = l:stepsize 
for j = l:stepsize 

L = i * Lmax/stepsize 
H = j * Hmax/stepsize 
objective(ij) = COMSOL calculations (L,H) 
end 

j =1 

end 

return objective 

This solution routine is shown graphically in Fig. 4. First, we 
calculate the extrema of the objectives that can be reached given 
a set of allowable domain dimensions. Function 1 (FI) passes this 
information into COMSOL where the temperature distribution for 
each case is calculated. As a result, we determine the work output, 
viscous resistance and the heat fluxes that correspond to each of 
the extreme domain sizes. This information is stored in a file. 

Next, we define a vector Vo that contains the initial guess for 
variables L and H. fminsearch passes this information into 
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Fig. 4. Flow diagram of the optimization routine, including the preprocessor. 


Function 2 (F2) which in turn calls COMSOL to perform temper¬ 
ature calculations for the specific domain. As a result, the “local” 
work output, resistance and heat fluxes are calculated. The values 
for the extrema calculated before are read from the file and used 
to normalize the objectives. If the algorithm determines that the 
current solution is optimal, it terminates the process, otherwise, 
it derives a new guess for domain dimensions which is passed 
back into the algorithm. 


5.2. Mesh dependence of the solution 

We conducted a sensitivity analysis in regard to the heat flux 
through the top surface of the domain (subject to convective heat 
flux). The mesh was refined in 5 steps, increasing the cell count by 
a factor of approximately 4 each time. The initial mesh contained 
516 cells. The largest mesh that COMSOL was able to solve con¬ 
tained over 132,000 cells. We found that relative to the total value 
of the heat flux, the error resulting from the coarse mesh is 
approximately 2%. Fig. 5 shows a plot of the heat flux post¬ 
processing data versus used mesh size. It shows that the step 
improvement is leveling off, thus a further increase of error can be 
assumed to be negligible. 


for work output, resistance, and heat fluxes in the feasible domain. 
Their values are dependent on the domain dimensions, where all 
other previously mentioned parameters (such as channel size, and 
temperature gradient) are assumed to be constant. This allows for 
visualization of the results using surface and contour plots. The 
calculation of the normalized heat fluxes requires COMSOL to be 
executed twice for each combination of domain width and height. 
This is very resource intensive, 1 and can severely slow down the 
entire process, especially when considering that the optimization is 
for only two decision variables. The optimization (the right hand 
side of Fig. 4) only requires one COMSOL evaluation. The hardware 
used is a computer equipped with a 3.2 GHz Pentium 4 processor, 
using 2 GB of RAM. This sensitivity analysis was done with grids of 
5 x 5, 10 x 10, 20 x 20, and 40 x 40 cells. The domain sizes were 
kept constant. The results are shown in Table 1. 

The sensitivity analysis in regard to the pre-processing and 
determination of the respective extreme values of the heat fluxes 
calculated by COMSOL shows that there is very little variation when 
the discretization of the grid is changed. This is because the func¬ 
tions for the heat fluxes are monotonic, and their maxima and 
minima occur on the boundary where L and H are at extrema as 
well. This is an important discovery, and could not have been 
assumed before doing this evaluation. 


A discretization of 40 x 40 cells requires 1600 COSMOL evaluations. 
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Fig. 6. Objective Function as a function of domain dimensions for PMMA. 


6. Results 

The feasible domain for a thermoacoustic stack is defined as 
0.1-0.5 m. This range is large, but was chosen in order to 
illustrate the behavior of the objective functions outside of the 
“practical” range. The results are presented in the form of 
contour plots of the objective function and a variation of the 
stack material, a variation of the initial guess (for constant 
material and weights), and finally a variation of the weights for 
each of the four objectives. 

6.2. Shape of the objective function 

As a result of the preprocessor, we can determine the overall 
shape of the objective function. Since it is only dependent on two 
parameters, it can be illustrated using a surface plot. One of these 
plots is shown in Fig. 6. This figure represents the objectives 
weighted equally at 25% each. 

In order to provide more insight in the behavior of this 
function depending on the material used for the stack, we ran 
this optimization for five additional (feasible) materials. The 
physical properties of the materials used are given in Table 2. 



Fig. 7. Behavior of the Optimization as a function of initial guess. 



Fig. 8. Contours for Case 5: Emphasis on the Heat Flux Objectives (Convection and 
Conduction). 

These results are given in the Appendix in Fig. 9. The overall 
trend is the same regardless of material, but we can notice slight 
variation in the shape of the objective function. The highly 
thermally conductive materials (aluminum and copper) are 
shown to be more sensitive to a variation in stack dimension 
than the alternative materials. 

6.2. Optimization behavior 

The behavior of the optimization routine can be shown when we 
consider a contour plot of the objective function and plot the 
domain size for each function call. This result is shown in Fig. 7. This 
shows how the optimizer changes the domain variables for each 
iteration step for different initial guesses (constant L, H varying). 
Because of the shape of the objective function in the area of the 
initial guesses, the solution converges to the same combination of L 
and H for each case. Fig. 10 in the Appendix illustrates additional 
cases of this analysis. 

6.3. Varying the weights 

For the PMMA 2 case, we varied the weights, all cases high¬ 
lighting emphasis on one objective or a group of objectives (power 
and resistance is one group, the heat fluxes another group). Table 3 
shows the weight distribution for all cases. 

Fig. 8 shows the case for an emphasis on just the heat loss 
via convection and conduction (Case 5). Multiplied by the open 
area coefficient, as defined in Equation (3.10). We can see that in 
order to minimize the heat losses, we must design the stack to 
be as small as possible. In addition of only focussing on the heat 
losses, we also investigated other weight distributions. The 
corresponding contour plots are shown in Fig. 11 in the 
Appendix. For example, disregarding the thermal losses 
completely, the emphasis on an increased height (that is greater 
cross sectional area of the stack) is apparent. Emphasizing solely 
on power, the design is required to be as large as possible 
(keeping in mind that in this case, the objective is to minimize 
the negative value). 


2 The tendencies revealed for PMMA are very similar for all other materials 
considered. 
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Table 2 

Properties of the materials used in the material variation. 


Material 

Density [kg/m 3 ] 

Heat capacity Q/kgK] 

Thermal conductivity [W/mK] 

Aluminium 

2700 

900 

160 

Copper 

8700 

385 

400 

Titanium 

4940 

710 

7.5 

Steel 

7850 

475 

44.5 

PMMA 

1190 

1420 

0.19 

Si0 2 

220 

730 

1.4 


Table 3 

Weight distribution (in %) for all objectives, corresponding to Fig. 11. 



7. Discussion and conclusions 

We have developed a basic model of a stack as part of a ther¬ 
moacoustic engine. This model is built in order to develop an 
understanding of the importance of the thermal losses encountered 
in the operation of thermoacoustic engines. As a design supple¬ 
ment, we included the calculation of the heat fluxes into an opti¬ 
mization routine, a method that is currently underused in 
thermoacoustics. In order to illustrate a dependency on the stack 
bulk material, we provided results for several different feasible 
stack materials. 

7.2. Discussion 

The surface plots for all materials show that the objective 
function is not convex. For this reason, there are several local 
optima but from a visual inspection, only one global optimum. It is 
obvious that the total length is always decreased to its imposed 
minimum value. The penalty functions keep this boundary from 
being violated. This location corresponds to a mechanically feasible 
design as it shows that the radius of the stack should be roughly 
20% larger than the length of the stack (again, provided the 
designer weighs all objectives equally). As a result of the material 
variation, we can identify slight changes in the absolute values of 
the objective function, but not a drastic change in the general 
tendencies and shape. The ideal design remains at the smallest 
value for L and a slightly larger value for the radius. 

The illustration of the optimization behavior showed all 
combinations of L and H that fminsearch considered. It shows that 
for most starting points, the area of the minimum objective value 
was reached. For starting points that were in the plateau regions of 
the objective function, fminsearch was not able to reach this global 
optimum, which should be a result of limited step size. For large 
initial values of L and D, the optimization found the other 
“optimum” which would result in as large a design as possible. This 
is, again, rooted in the limited initial step size. 

The variation of the weights for each objective shows how the 
design would change if emphasis is given to one objective in 
particular, for example as much power while essentially ignoring 
all losses. In this case, the result showed that as large a design as 
possible is the solution (which we could have derived from inves¬ 
tigating the function for power in Equation (3.7). A more interesting 
derivation is the design with an emphasis on avoiding the 
convective/conductive heat loss or the conductive heat loss (and 


the combination of both) which can be interpreted as “as small as 
possible.” The following design criteria can then be deduced: 

• Design as small as possible to minimize thermal losses; 

• Increase overall stack radius to increase power output, at 
cost of higher thermal losses, with optimum at about 
H « 1.2L; 

• Design as large as possible to maximize power output. 

Obviously, choice 1 and choice 3 are conflicting. This illustrates 
the tradeoff that we have to consider when designing a thermoa¬ 
coustic engine. Ultimately, it is a personal choice whether to 
emphasize the thermal behavior or the acoustic behavior. In 
closing, we will now elaborate briefly on how to improve the 
current optimization scheme in order to provide more general 
design information. 

7.2. Conclusions 

The presented results are a good first approach to using optimi¬ 
zation principles in the design of thermoacoustic devices. Considering 
four objectives simultaneously has shown that for a given length, 
there is an optimal diameter that results in the least amount of heat 
loss across the boundary surface and through the cooling water (that 
is for the case of equal weights given to all four objectives). We have 
determined a design statement for a driving stack of a thermoacoustic 
engine. In the future, this work has to be expanded to include ther¬ 
moacoustic refrigerators, as this technology is the primary application 
of thermoacoustics with the most potential impact. 

In order to fully represent a thermoacoustic system, we have to 
consider the phasing relationship between pressure and velocity and 
its impact on heat transfer. This would result in an added influence of 
the channel size, which is currently assumed to be a constant. The 
effect of delayed heat transfer, which, with regard to standing wave 
engines, is an important influence, is neglected. In standing wave 
engines, the heat transfer from gas to solid has to be delayed in order 
to achieve the correct phasing between pressure and velocity oscil¬ 
lations. This delay is shorter the smaller the channel is. If the channel 
size is equal to or less than the thermal penetration depth, there is 
ideal contact, i.e. no delay in the heat transfer. Ultimately, a consid¬ 
eration of acoustic theory and improved models of thermoacoustic 
work and loss mechanisms have to be included. 

Acknowledgements 

This work was supported through the NSF East Asia and Pacific 
Summer Institute 2008 (OISE-0812781). 

Appendix 

Plots corresponding to the material variation 

The following surface plots illustrate the effect of the variation of 
the material used for the calculations in COMSOL. This is a combi¬ 
nation of all four objectives, with equal weights of 25%. For the 
materials with smaller thermal conductivity, the objective function 
is more flat than for the metal cases. 

Plots showing the behavior of ‘fminsearch ” 

This plot illustrates the behavior of the optimizer for different 
starting points. Each of the five cases shows different starting 
points for a constant domain length L. The material used in PMMA, 
and the objective function is shown as a contour plot. 
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Fig. 9. Variation of the stack material. 
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H 

Varying H with L = 0.5 


Fig. 10. Variation of the starting points for constant L each (Material PMMA). 
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H H 

Case 1: Emphasis on Power Case 2: Emphasis on Resistance 



H H 

Case 3: Emphasis on Convective Heat Flux Case 4: Emphasis on Conductive Heat Flux 



H H 

Case 5: Emphasis on Heat Loss Case 6: Emphasis on Power and Resistance 


Fig. 11. Variation of the weights for each objective. 
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